%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%simulate model thru epidemic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[neglogpost,cmat,comat,cymat,cons_monthly,pidy_vec_beliefs,pido_vec_beliefs,par_pi1,par_pi2,bprimemat,Umat]=simulate_model(guess,Ivec,muvec,par,U,Y,mpar,grid,solflag,Umat,bprimemat,Sy,Ry,Iy,So,Ro,Io,newdeaths_weekly,dat_young,dat_old,dt_cons,fminbnd_options,c_NoEpiSim)


par.muvec_scale=abs(guess(1));
par.pidy_ini=abs(guess(2));
par.pido_ini=abs(guess(3));
par.gyoung=abs(guess(4));
par.gold=abs(guess(5));
par.kappa=abs(guess(6));


try
    %back out pi1 and pi2 from R0 and kappa
    %kappa=(par.pi1*(par.cry_noepi_ini*par.sy_shr+par.cro_noepi_ini*(1-par.sy_shr)))/(par.pi1*(par.cry_noepi_ini*par.sy_shr+par.cro_noepi_ini*(1-par.sy_shr))+par.pi2); %take weighted avg of old and young consumption
    %par.Rnot=(par.pi2+par.pi1*par.cry_noepi_ini*par.sy_shr+par.pi1*par.cro_noepi_ini*(1-par.sy_shr) ) /(par.piry*par.sy_shr+par.piro*(1-par.sy_shr)+par.pidy*par.sy_shr+par.pido*(1-par.sy_shr));
    tmp=(par.cry_noepi_ini*par.sy_shr+par.cro_noepi_ini*(1-par.sy_shr));
    tmp2=(par.piry*par.sy_shr+par.piro*(1-par.sy_shr)+par.pidy*par.sy_shr+par.pido*(1-par.sy_shr));
    par.pi2=(1-par.kappa)*par.R0*tmp2;
    par.pi1=par.kappa*par.pi2/((1-par.kappa)*tmp);
    par_pi1=par.pi1;
    par_pi2=par.pi2;
    
    %scale containment measure
    muvec=muvec*par.muvec_scale;
    
    %constant gain belief evolution
    par.pidy_vec_beliefs=Ivec*NaN; %initial cfr
    par.pido_vec_beliefs=Ivec*NaN; %initial cfr
    
    par.pidy_vec_beliefs(1)=par.pidy_ini;
    par.pido_vec_beliefs(1)=par.pido_ini;
    
    for rr=2:1:numel(Ivec)
        par.pidy_vec_beliefs(rr)=par.pidy_vec_beliefs(rr-1)+par.gyoung*(par.pidy_vec_true(rr)-par.pidy_vec_beliefs(rr-1));
        par.pido_vec_beliefs(rr)=par.pido_vec_beliefs(rr-1)+par.gold*(par.pido_vec_true(rr)-par.pido_vec_beliefs(rr-1));
    end
    
    %time varying beliefs of probabilities of dying and recovering
    par.pidy_vec_beliefs=par.pidy_vec_beliefs;%*0+par.pidy_vec_true;%*0+7*0.001/par.days;
    par.piry_vec_beliefs=7*1/par.days-par.pidy_vec_beliefs;
    par.pido_vec_beliefs=par.pido_vec_beliefs;%*0+par.pido_vec_true;%*0+7*0.035/par.days;
    par.piro_vec_beliefs=7*1/par.days-par.pido_vec_beliefs;
    
    if par.piry_vec_beliefs<0,error('piry <0');end
    if par.piro_vec_beliefs<0,error('piro <0');end
    
    
    %Store original variables
    Ivec_base=Ivec;
    muvec_base=muvec;
    
    par.pidy_vec_beliefs_base=par.pidy_vec_beliefs;
    par.piry_vec_beliefs_base=par.piry_vec_beliefs;
    par.pido_vec_beliefs_base=par.pido_vec_beliefs;
    par.piro_vec_beliefs_base=par.piro_vec_beliefs;
    
    Umat_base=Umat;
    bprimemat_base=bprimemat;
    
    
    if mpar.infect_surprise==0 %perfect foresight for infections and containment
        
        %simulate value functions of agents by age and health types
        sim_val_fun;
        
        %store assets first wave
        bprimemat_foresight=bprimemat;
        
    elseif mpar.infect_surprise==1  %simulate first wave until period 21 (assume steady state in period 22). Then surprise in period 22 that second wave takes place.
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%First wave
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %first wave -- take relevant stuff for first wave
        Ivec=Ivec(1:mpar.infect_surprise_week);
        Ivec=[Ivec,0];
        
        muvec=muvec(1:mpar.infect_surprise_week+1);
        
        par.pidy_vec_beliefs=par.pidy_vec_beliefs(1:mpar.infect_surprise_week+1);
        par.piry_vec_beliefs=par.piry_vec_beliefs(1:mpar.infect_surprise_week+1);
        par.pido_vec_beliefs=par.pido_vec_beliefs(1:mpar.infect_surprise_week+1);
        par.piro_vec_beliefs=par.piro_vec_beliefs(1:mpar.infect_surprise_week+1);
        
        Umat=Umat_base;
        bprimemat=bprimemat_base;
        
        Umat.ro=Umat.ro(:,end-mpar.infect_surprise_week:end);
        Umat.ry=Umat.ry(:,end-mpar.infect_surprise_week:end);
        Umat.io=Umat.io(:,end-mpar.infect_surprise_week:end);
        Umat.iy=Umat.iy(:,end-mpar.infect_surprise_week:end);
        Umat.so=Umat.so(:,end-mpar.infect_surprise_week:end);
        Umat.sy=Umat.sy(:,end-mpar.infect_surprise_week:end);
        
        bprimemat.ro=bprimemat.ro(:,end-mpar.infect_surprise_week:end);
        bprimemat.ry=bprimemat.ry(:,end-mpar.infect_surprise_week:end);
        bprimemat.io=bprimemat.io(:,end-mpar.infect_surprise_week:end);
        bprimemat.iy=bprimemat.iy(:,end-mpar.infect_surprise_week:end);
        bprimemat.so=bprimemat.so(:,end-mpar.infect_surprise_week:end);
        bprimemat.sy=bprimemat.sy(:,end-mpar.infect_surprise_week:end);
        
        %simulate value functions of agents by age and health types
        sim_val_fun;
        
        %store assets and Utility first wave
        bprimemat_first_wave=bprimemat;
        Umat_first_wave=Umat;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Second wave
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %re-set variables to base
        Ivec=Ivec_base;
        muvec=muvec_base;
        
        par.pidy_vec_beliefs=par.pidy_vec_beliefs_base;
        par.piry_vec_beliefs=par.piry_vec_beliefs_base;
        par.pido_vec_beliefs=par.pido_vec_beliefs_base;
        par.piro_vec_beliefs=par.piro_vec_beliefs_base;
        
        Umat=Umat_base;
        bprimemat=bprimemat_base;
        
        %second wave -- take relevant stuff for second wave
        Ivec=Ivec(mpar.infect_surprise_week+1:end);
        muvec=muvec(mpar.infect_surprise_week+1:end);
        
        par.pidy_vec_beliefs=par.pidy_vec_beliefs(mpar.infect_surprise_week+1:end);
        par.piry_vec_beliefs=par.piry_vec_beliefs(mpar.infect_surprise_week+1:end);
        par.pido_vec_beliefs=par.pido_vec_beliefs(mpar.infect_surprise_week+1:end);
        par.piro_vec_beliefs=par.piro_vec_beliefs(mpar.infect_surprise_week+1:end);
        
        Umat.ro=Umat.ro(:,mpar.infect_surprise_week+1:end);
        Umat.ry=Umat.ry(:,mpar.infect_surprise_week+1:end);
        Umat.io=Umat.io(:,mpar.infect_surprise_week+1:end);
        Umat.iy=Umat.iy(:,mpar.infect_surprise_week+1:end);
        Umat.so=Umat.so(:,mpar.infect_surprise_week+1:end);
        Umat.sy=Umat.sy(:,mpar.infect_surprise_week+1:end);
        
        bprimemat.ro=bprimemat.ro(:,mpar.infect_surprise_week+1:end);
        bprimemat.ry=bprimemat.ry(:,mpar.infect_surprise_week+1:end);
        bprimemat.io=bprimemat.io(:,mpar.infect_surprise_week+1:end);
        bprimemat.iy=bprimemat.iy(:,mpar.infect_surprise_week+1:end);
        bprimemat.so=bprimemat.so(:,mpar.infect_surprise_week+1:end);
        bprimemat.sy=bprimemat.sy(:,mpar.infect_surprise_week+1:end);
        
        %simulate value functions of agents by age and health types
        sim_val_fun;
        
        %re-set variables to base
        Ivec=Ivec_base;
        muvec=muvec_base;
        
        par.pidy_vec_beliefs=par.pidy_vec_beliefs_base;
        par.piry_vec_beliefs=par.piry_vec_beliefs_base;
        par.pido_vec_beliefs=par.pido_vec_beliefs_base;
        par.piro_vec_beliefs=par.piro_vec_beliefs_base;
        
        %store assets and utility second wave
        bprimemat_second_wave=bprimemat;
        Umat_second_wave=Umat;
        
        %put together assets and utility from both waves
        bprimemat.ro=[bprimemat_first_wave.ro(:,1:mpar.infect_surprise_week) bprimemat_second_wave.ro];
        bprimemat.ry=[bprimemat_first_wave.ry(:,1:mpar.infect_surprise_week) bprimemat_second_wave.ry];
        bprimemat.io=[bprimemat_first_wave.io(:,1:mpar.infect_surprise_week) bprimemat_second_wave.io];
        bprimemat.iy=[bprimemat_first_wave.iy(:,1:mpar.infect_surprise_week) bprimemat_second_wave.iy];
        bprimemat.so=[bprimemat_first_wave.so(:,1:mpar.infect_surprise_week) bprimemat_second_wave.so];
        bprimemat.sy=[bprimemat_first_wave.sy(:,1:mpar.infect_surprise_week) bprimemat_second_wave.sy];
        
        Umat.ro=[Umat_first_wave.ro(:,1:mpar.infect_surprise_week) Umat_second_wave.ro];
        Umat.ry=[Umat_first_wave.ry(:,1:mpar.infect_surprise_week) Umat_second_wave.ry];
        Umat.io=[Umat_first_wave.io(:,1:mpar.infect_surprise_week) Umat_second_wave.io];
        Umat.iy=[Umat_first_wave.iy(:,1:mpar.infect_surprise_week) Umat_second_wave.iy];
        Umat.so=[Umat_first_wave.so(:,1:mpar.infect_surprise_week) Umat_second_wave.so];
        Umat.sy=[Umat_first_wave.sy(:,1:mpar.infect_surprise_week) Umat_second_wave.sy];
        
    end
    
    %get consumption policy functions
    cmat.ro=(par.w+(1+par.r)*grid.b(mpar.grid_sim_idx)'-bprimemat.ro);
    cmat.ry=(par.w+(1+par.r)*grid.b(mpar.grid_sim_idx)'-bprimemat.ry);
    cmat.io=(par.w+(1+par.r)*grid.b(mpar.grid_sim_idx)'-bprimemat.io);
    cmat.iy=(par.w+(1+par.r)*grid.b(mpar.grid_sim_idx)'-bprimemat.iy);
    cmat.so=(par.w+(1+par.r)*grid.b(mpar.grid_sim_idx)'-bprimemat.so);
    cmat.sy=(par.w+(1+par.r)*grid.b(mpar.grid_sim_idx)'-bprimemat.sy);
    
    
    %simulate consumption, assets etc over time starting in epi t=0 for
    %some given initial assets (i.e. track young and old at initial assets
    %through epi)
    
    ZZ=numel(Umat.so(1,:)); %periods to simulate
    
    c_EpiSim.so=NaN*zeros(ZZ,1);c_EpiSim.sy=NaN*zeros(ZZ,1);
    c_EpiSim.ro=NaN*zeros(ZZ,1);c_EpiSim.ry=NaN*zeros(ZZ,1);
    c_EpiSim.io=NaN*zeros(ZZ,1);c_EpiSim.iy=NaN*zeros(ZZ,1);
    
    b_EpiSim.sy=NaN*zeros(ZZ,1);b_EpiSim.so=NaN*zeros(ZZ,1);
    b_EpiSim.iy=NaN*zeros(ZZ,1);b_EpiSim.io=NaN*zeros(ZZ,1);
    b_EpiSim.ry=NaN*zeros(ZZ,1);b_EpiSim.ro=NaN*zeros(ZZ,1);
    
    b_EpiSim.so(1)=grid.b(par.grid_b_old_assets_idx);%initial assets, old
    b_EpiSim.sy(1)=grid.b(par.grid_b_young_assets_idx);%initial assets, young
    b_EpiSim.ro(1)=grid.b(par.grid_b_old_assets_idx);%initial assets, old
    b_EpiSim.ry(1)=grid.b(par.grid_b_young_assets_idx);%initial assets, young
    b_EpiSim.io(1)=grid.b(par.grid_b_old_assets_idx);%initial assets, old
    b_EpiSim.iy(1)=grid.b(par.grid_b_young_assets_idx);%initial assets, young
    
    %notation:
    %grid.b: grid of assets, t-1
    %bprime: grid of assets, t (conditional on t-1 assets), policy function
    %for assets
    %c: policy function for consumption (budget constraint) as a function
    %of assets, t-1 and t
    
    
    for jj=2:1:ZZ+1
        b_EpiSim.so(jj) = interp1(grid.b(mpar.grid_sim_idx),bprimemat.so(:,jj-1),b_EpiSim.so(jj-1),'linear','extrap'); %get path of assets given initial assets
        c_EpiSim.so(jj-1) = (par.w+(1+par.r)*b_EpiSim.so(jj-1)-b_EpiSim.so(jj)); %use budget to calc. cons.
        
        b_EpiSim.sy(jj) = interp1(grid.b(mpar.grid_sim_idx),bprimemat.sy(:,jj-1),b_EpiSim.sy(jj-1),'linear','extrap'); %get path of assets given initial assets
        c_EpiSim.sy(jj-1) = (par.w+(1+par.r)*b_EpiSim.sy(jj-1)-b_EpiSim.sy(jj)); %use budget to calc. cons.
        
        b_EpiSim.ro(jj) = interp1(grid.b(mpar.grid_sim_idx),bprimemat.ro(:,jj-1),b_EpiSim.ro(jj-1),'linear','extrap'); %get path of assets given initial assets
        c_EpiSim.ro(jj-1) = (par.w+(1+par.r)*b_EpiSim.ro(jj-1)-b_EpiSim.ro(jj)); %use budget to calc. cons.
        
        b_EpiSim.ry(jj) = interp1(grid.b(mpar.grid_sim_idx),bprimemat.ry(:,jj-1),b_EpiSim.ry(jj-1),'linear','extrap'); %get path of assets given initial assets
        c_EpiSim.ry(jj-1) = (par.w+(1+par.r)*b_EpiSim.ry(jj-1)-b_EpiSim.ry(jj)); %use budget to calc. cons.
        
        b_EpiSim.io(jj) = interp1(grid.b(mpar.grid_sim_idx),bprimemat.io(:,jj-1),b_EpiSim.io(jj-1),'linear','extrap'); %get path of assets given initial assets
        c_EpiSim.io(jj-1) = (par.w+(1+par.r)*b_EpiSim.io(jj-1)-b_EpiSim.io(jj)); %use budget to calc. cons.
        
        b_EpiSim.iy(jj) = interp1(grid.b(mpar.grid_sim_idx),bprimemat.iy(:,jj-1),b_EpiSim.iy(jj-1),'linear','extrap'); %get path of assets given initial assets
        c_EpiSim.iy(jj-1) = (par.w+(1+par.r)*b_EpiSim.iy(jj-1)-b_EpiSim.iy(jj)); %use budget to calc. cons.
    end    
    
    %per capita consumption of young and old starting epi with initial
    %assets
    cymat_new=(c_EpiSim.sy'.*Sy'+c_EpiSim.iy'.*Iy'+c_EpiSim.ry'.*Ry')./(Sy'+Iy'+Ry');
    comat_new=(c_EpiSim.so'.*So'+c_EpiSim.io'.*Io'+c_EpiSim.ro'.*Ro')./(So'+Io'+Ro');
    
    %padding no epi consumption path for the periods before the start
    %of the simulation
    cymat=[c_NoEpiSim.sy(1:mpar.timestart-1)' cymat_new];%note that c_NoEpiSim.sy is per capita and that I=R=0, S=1.
    comat=[c_NoEpiSim.so(1:mpar.timestart-1)' comat_new]; %note that c_NoEpiSim.so is per capita and that I=R=0, S=1.
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Epi Simulation results: consumption old and young
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %express epi consumption as percent of no epi consumption baseline
    cytmp=(cymat./c_NoEpiSim.sy(1:numel(cymat))'-1)*100;%note that c_NoEpiSim.sy is per capita and that I=R=0, S=1.
    cotmp=(comat./c_NoEpiSim.so(1:numel(cymat))'-1)*100;%note that c_NoEpiSim.so is per capita and that I=R=0, S=1.
    
    % Put data into timetable
    cons_weekly.Cy=[cytmp(1:numel(Ivec)+mpar.timestart-1)'];
    cons_weekly.Co=[cotmp(1:numel(Ivec)+mpar.timestart-1)'];
    
    tty = timetable(newdeaths_weekly.Time,cons_weekly.Cy,'VariableNames',{'Consy'});
    tto = timetable(newdeaths_weekly.Time,cons_weekly.Co,'VariableNames',{'Conso'});
    
    cons_monthly.Cy=retime(tty,'monthly','mean');
    cons_monthly.Co=retime(tto,'monthly','mean');
    
    %March 1 2020 thru April 30 2021
    disty=(dat_young(1:14)-cons_monthly.Cy.Consy(1:14)')/100; %inverse V weighting matrix is in original units (i.e. not in percent)
    disto=(dat_old(1:14)-cons_monthly.Co.Conso(1:14)')/100;
    
    dist=[disty';disto'];
    criterion=dist'*par.invVVmat*dist;
    loglik=numel(dist)/2*log(1/2/pi)-1/2*par.logdetVVmat-1/2*criterion; %in logs
    
    %uniform priors for all parameters
    logprior=0;
    
    %par.guess_name=[{'mu_scale   '};{'pid0_young'};{'pid0_old  '};{'g_young    '};{'g_old      '};{'pi1shr     '};{'R0ini     '}];%{'pi1shr     '};{'R0ini     '}]; %];
    for gg=1:1:numel(guess)
        if gg==1 %muscale -- uniform prior
            logprior=logprior+log(unifpdf(abs(guess(gg)),par.LB(gg),par.UB(gg)));
            
        elseif gg==2 %pid0_young -- uniform prior
            logprior=logprior+log(unifpdf(abs(guess(gg)),par.LB(gg),par.UB(gg)));
            
        elseif gg==3 %pid0_old -- uniform prior
            logprior=logprior+log(unifpdf(abs(guess(gg)),par.LB(gg),par.UB(gg)));
            
        elseif gg==4 %g_young -- uniform prior
            logprior=logprior+log(unifpdf(abs(guess(gg)),par.LB(gg),par.UB(gg)));
            
        elseif gg==5 %g_old -- uniform prior
            logprior=logprior+log(unifpdf(abs(guess(gg)),par.LB(gg),par.UB(gg)));
            
        elseif gg==6 %kappa -- uniform prior
            logprior=logprior+log(unifpdf(abs(guess(gg)),par.LB(gg),par.UB(gg)));
        end
    end
    
    logpost=loglik+logprior;
    
    neglogpost=-logpost; %flip sign for maximization using fmincon
    
    %if isinf(logpost)
    %    keyboard
    %end
    
    %for output of this function
    pidy_vec_beliefs=par.pidy_vec_beliefs;
    pido_vec_beliefs=par.pido_vec_beliefs;
    
catch
    
    neglogpost=1e5;
    %keyboard
    
end